********************************************************************************
********************************************************************************
** Robustness analysis 3
** Potentially confounding factors when switching type of school.
** Includes changes in:
**	Private tuition.
**	Child labour.
**	Child health.
**	General life satisfaction in the household.
**	Household consumption.
********************************************************************************
********************************************************************************

preserve

** Merge variables together

* Load education module from Y3
use "$dataraw_path\nps\2012-2013\HH_SEC_C.dta", clear

keep y3_hhid indidy3 hh_c05 hh_c09 hh_c12 hh_c28_1-hh_c28_8

merge 1:m y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\NPSY3.PANEL.KEY.dta", ///
	keepusing(y2_hhid indidy2)
drop if _merge==2
drop if indidy2==.
drop _merge

gen enrolled_y3 = (hh_c05==1) if hh_c05!=.
keep if enrolled_y3==1
gen grade_y3 = hh_c09
gen private_y3 = (!inlist(hh_c12,1,2)) if hh_c12!=.
drop if private_y3==.
gen exp1_y3 = hh_c28_1
gen exp2_y3 = hh_c28_2
gen exp3_y3 = hh_c28_3
gen exp4_y3 = hh_c28_4
gen exp5_y3 = hh_c28_5
gen exp6_y3 = hh_c28_6
gen exp7_y3 = hh_c28_7
gen exp8_y3 = hh_c28_8
drop hh_c05 hh_c09 hh_c12 hh_c28_1-hh_c28_8 enrolled_y3

* Merge with education module from Y2
merge 1:1 y2_hhid indidy2 using "$dataraw_path\nps\2010-2011\HH_SEC_C", ///
	keepusing(hh_c05 hh_c09 hh_c12 hh_c28_1-hh_c28_8)
drop if _merge==2
drop _merge

gen enrolled_y2 = (hh_c05==1) if hh_c05!=.
drop if enrolled_y2==0 | enrolled_y2==.
gen grade_y2 = hh_c09
gen private_y2 = (!inlist(hh_c12,1,2)) if hh_c12!=.
drop if private_y2==.
gen exp1_y2 = hh_c28_1
gen exp2_y2 = hh_c28_2
gen exp3_y2 = hh_c28_3
gen exp4_y2 = hh_c28_4
gen exp5_y2 = hh_c28_5
gen exp6_y2 = hh_c28_6
gen exp7_y2 = hh_c28_7
gen exp8_y2 = hh_c28_8
drop hh_c05 hh_c09 hh_c12 hh_c28_1-hh_c28_8 enrolled_y2

* Keep only students attending primary or lower/upper secondary school
drop if inlist(grade_y2,1,2,19,20) | grade_y2>=25
drop if inlist(grade_y3,1,2,19,20) | grade_y3>=25

tab private_y2 private_y3, row

* Merge with child health and average household subjective welfare from Y3
merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_D", ///
	keepusing(hh_d10 hh_d14)
drop if _merge==2
drop _merge

gen hospitalization_y3 = (hh_d10==1 | hh_d14==1) if hh_d10!=. | hh_d14!=.

merge m:m y3_hhid using "$dataraw_path\nps\2012-2013\HH_SEC_G", ///
	keepusing(hh_g03_1 hh_g03_6 hh_g03_8)
drop if _merge==2
drop _merge

replace hh_g03_1 = . if hh_g03_1==8
replace hh_g03_6 = . if hh_g03_6==8
replace hh_g03_8 = . if hh_g03_8==8

gen sat_temp1 = (inlist(hh_g03_1,1,2,3)) if hh_g03_1!=.
gen sat_temp2 = (inlist(hh_g03_6,1,2,3)) if hh_g03_6!=.
gen sat_temp3 = (inlist(hh_g03_8,1,2,3)) if hh_g03_8!=.

bys y3_hhid: egen sat_health1_y3 = mean(hh_g03_1)
bys y3_hhid: egen sat_health2_y3 = mean(sat_temp1)
bys y3_hhid: egen sat_educ1_y3 = mean(hh_g03_6)
bys y3_hhid: egen sat_educ2_y3 = mean(sat_temp2)
bys y3_hhid: egen sat_life1_y3 = mean(hh_g03_8)
bys y3_hhid: egen sat_life2_y3 = mean(sat_temp3)

duplicates drop y3_hhid indidy3, force
drop hh_d10 hh_d14 hh_g03_1 hh_g03_6 hh_g03_8 sat_temp1 sat_temp2 sat_temp3

* Merge with child health and average household subjective welfare from Y2
merge 1:1 y2_hhid indidy2 using "$dataraw_path\nps\2010-2011\HH_SEC_D", ///
	keepusing(hh_d10 hh_d14)
drop if _merge==2
drop _merge

gen hospitalization_y2 = (hh_d10==1 | hh_d14==1) if hh_d10!=. | hh_d14!=.

merge m:m y2_hhid using "$dataraw_path\nps\2010-2011\HH_SEC_G", ///
	keepusing(hh_g03_1 hh_g03_6 hh_g03_8)
drop if _merge==2
drop _merge

replace hh_g03_1 = . if hh_g03_1==8
replace hh_g03_6 = . if hh_g03_6==8
replace hh_g03_8 = . if hh_g03_8==8

gen sat_temp1 = (inlist(hh_g03_1,1,2,3)) if hh_g03_1!=.
gen sat_temp2 = (inlist(hh_g03_6,1,2,3)) if hh_g03_6!=.
gen sat_temp3 = (inlist(hh_g03_8,1,2,3)) if hh_g03_8!=.

bys y2_hhid: egen sat_health1_y2 = mean(hh_g03_1)
bys y2_hhid: egen sat_health2_y2 = mean(sat_temp1)
bys y2_hhid: egen sat_educ1_y2 = mean(hh_g03_6)
bys y2_hhid: egen sat_educ2_y2 = mean(sat_temp2)
bys y2_hhid: egen sat_life1_y2 = mean(hh_g03_8)
bys y2_hhid: egen sat_life2_y2 = mean(sat_temp3)

duplicates drop y2_hhid indidy2, force
drop hh_d10 hh_d14 hh_g03_1 hh_g03_6 hh_g03_8 sat_temp1 sat_temp2 sat_temp3

* Child work and parental work from Y3 (cannot include self-employed hours as it is not reported in Y2)
merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_E", ///
	keepusing(hh_e08c hh_e32 hh_e50 hh_e64 hh_e69)
drop if _merge==2
drop _merge

egen work_child1_y3 = rowtotal(hh_e32 hh_e50 hh_e64 hh_e69) // student work hours (intensive margin)
gen work_child2_y3 = (work_child1_y3>=1 | hh_e08c==1) // student work dummy (extensive margin)

drop hh_e08c-hh_e69

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_B", ///
	keepusing(hh_b04) nogen // get age of all household members
rename hh_b04 age

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_C", ///
	keepusing(hh_c05) nogen // get school enrolment for all household members

merge 1:1 y3_hhid indidy3 using "$dataraw_path\nps\2012-2013\HH_SEC_E", ///
	keepusing(hh_e08c hh_e32 hh_e50 hh_e64 hh_e69) nogen // get labor information for all household members

egen work_temp1 = rowtotal(hh_e32 hh_e50 hh_e64 hh_e69) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work hours 
replace work_temp1 = 168 if work_temp1!=. & work_temp1>168
gen work_temp2 = (work_temp1>=1 & work_temp1!=. | hh_e08c==1) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work dummy

bys y3_hhid: egen work_all1_y3 = mean(work_temp1) // average work hours in household
bys y3_hhid: egen work_all2_y3 = mean(work_temp2) // share working in household

drop hh_c05-work_temp2 age
drop if y2_hhid==""


* Child work and parental work from Y2 (cannot include self-employed hours as it is not reported in Y2)
merge 1:1 y2_hhid indidy2 using "$dataraw_path\nps\2010-2011\HH_SEC_E1", ///
	keepusing(hh_e51 hh_e25 hh_e40 hh_e75 hh_e78)
drop if _merge==2
drop _merge

egen work_child1_y2 = rowtotal(hh_e25 hh_e40 hh_e75 hh_e78) // student work hours (intensive margin)
gen work_child2_y2 = (work_child1_y2>=1 | hh_e51==1) // student work dummy (extensive margin)

drop hh_e25-hh_e78

merge 1:1 y2_hhid indidy2 using "$dataraw_path\nps\2010-2011\HH_SEC_B", ///
	keepusing(hh_b04) nogen // get age of all household members
rename hh_b04 age

merge 1:1 y2_hhid indidy2 using "$dataraw_path\nps\2010-2011\HH_SEC_C", ///
	keepusing(hh_c05) nogen // get school enrolment for all household members

merge 1:1 y2_hhid indidy2 using "$dataraw_path\nps\2010-2011\HH_SEC_E1", ///
	keepusing(hh_e51 hh_e25 hh_e40 hh_e75 hh_e78) nogen // get labor information for all household members

egen work_temp1 = rowtotal(hh_e25 hh_e40 hh_e75 hh_e78) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work hours 
replace work_temp1 = 168 if work_temp1!=. & work_temp1>168
gen work_temp2 = (work_temp1>=1 & work_temp1!=. | hh_e51==1) ///
	if age>=15 & age!=. & hh_c05!=1 // non-student household members work dummy

bys y2_hhid: egen work_all1_y2 = mean(work_temp1) // average work hours in household
bys y2_hhid: egen work_all2_y2 = mean(work_temp2) // share working in household

drop hh_c05-work_temp2 age
drop if y3_hhid==""


* Household consumption from Y3
merge m:1 y3_hhid using "$dataraw_path\nps\2012-2013\ConsumptionNPS3", ///
	keepusing(expmR adulteq)
drop if _merge==2
drop _merge

gen cons_y3 = expmR/adulteq

drop expmR adulteq


* Household consumption from Y2
merge m:1 y2_hhid using "$dataraw_path\nps\2010-2011\TZY2.HH.Consumption.dta", ///
	keepusing(expmR adulteq)
drop if _merge==2
drop _merge

gen cons_y2 = expmR/adulteq

drop expmR adulteq


** Create variables used in regressions

* Reshape from wide to long
reshape long grade_y private_y exp1_y exp2_y exp3_y exp4_y exp5_y exp6_y ///
exp7_y exp8_y hospitalization_y sat_health1_y sat_health2_y sat_educ1_y ///
sat_educ2_y sat_life1_y sat_life2_y work_child1_y work_child2_y work_all1_y ///
work_all2_y cons_y, i(y3_hhid-indidy2) j(round)

foreach var of varlist *_y {
   	local newname = substr("`var'", 1, length("`var'")-2)
   	rename `var' `newname'
}

* Create unique student ID
egen id = group(y3_hhid indidy3)

* Set panel and time variables
xtset id round

* Log of variables
gen exp5_log = log(1+exp5)
gen work_child1_log = log(1+work_child1)
gen work_all1_log = log(1+work_all1)
gen cons_log = log(1+cons)

* Sample dummy and empty matrix for sample means and standard deviations
gen sample_temp = 1 if exp5_log!=. & l.exp5_log!=. & sat_life1!=. & ///
	l.sat_life1!=. & cons_log!=.
bys id: egen sample = max(sample_temp)
sort id round

matrix a = J(2,6,.)


** Regress changes in outcome variables on changes in private status

* Money spent on private tuition (log)
eststo rob3_1: reg d.exp5_log d.private if sample==1, vce(robust) // log of tuition expenditures

sum exp5_log if sample==1 | f.sample==1
mat a[1,1] = `r(mean)'
mat a[2,1] = `r(sd)'

* Child work
eststo rob3_2: reg d.work_child2 d.private if sample==1, vce(robust) // dummy indicates whether the student worked
eststo rob3_3: reg d.work_child1_log d.private if sample==1, vce(robust) // log of hours worked by the student

sum work_child2 if sample==1 | f.sample==1
mat a[1,2] = `r(mean)'
mat a[2,2] = `r(sd)'

sum work_child1_log if sample==1 | f.sample==1
mat a[1,3] = `r(mean)'
mat a[2,3] = `r(sd)'

* Child health
eststo rob3_4: reg d.hospitalization d.private if sample==1, vce(robust) // dummy indicating whether the student has been hospitalized in past 12 months

sum hospitalization if sample==1 | f.sample==1
mat a[1,4] = `r(mean)'
mat a[2,4] = `r(sd)'


* Household subjective well-being
//reg d.sat_health1 d.private, vce(robust) // average household satisfaction with own health
//reg d.sat_health2 d.private, vce(robust) // share of household members satisfied with own health
//reg d.sat_educ1 d.private, vce(robust) // average household satisfaction with education available
//reg d.sat_educ2 d.private, vce(robust) // share of household members satisfied with education available
eststo rob3_5: reg d.sat_life1 d.private if sample==1, vce(robust) // average household satisfaction with own life
//reg d.sat_life2 d.private if sample==1, vce(robust) // share of household members satisfied with own life

sum sat_life1 if sample==1 | f.sample==1
mat a[1,5] = `r(mean)'
mat a[2,5] = `r(sd)'


* Consumption
eststo rob3_6: reg d.cons_log d.private if sample==1, vce(robust) // log of household consumption per adult-equivalent member

sum cons_log if sample==1 | f.sample==1
mat a[1,6] = `r(mean)'
mat a[2,6] = `r(sd)'


** Output
esttab rob3_* using "$out_path\tablea9.tex", ///
replace se stats(N r2, fmt(%12.3gc) labels("\(N\)" "\(R^2\)")) compress ///
nomtitles starlevels(* 0.1 ** 0.05 *** 0.01) substitute(\_ _) b(3) ///
k(D.private) varl(D.private "\$\Delta Private$")

file open rob3 using "$out_path\tablea9_stats.tex", write replace
file write rob3 "\textit{Sample mean} & " (round(a[1,1],0.001)) " & " ///
	(round(a[1,2],0.001)) " & " (round(a[1,3],0.001)) " & " ///
	(round(a[1,4],0.001)) " & " (round(a[1,5],0.001)) " & " /// 
	(round(a[1,6],0.001)) " \\ " _n ///
	"\textit{Sample std. dev.} & " (round(a[2,1],0.001)) " & " ///
	(round(a[2,2],0.001)) " & " (round(a[2,3],0.001)) " & " ///
	(round(a[2,4],0.001)) " & " (round(a[2,5],0.001)) " & " ///
	(round(a[2,6],0.001)) " \\"
file close rob3

restore


